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ABSTRACT 

Shear modulus of solid neutron star crust is calculated by tliermodynamic perturbation 
theory taking into account ion motion. At given density the crust is modelled as a 
body-centered cubic Coulomb crystal of fully ionized atomic nuclei of one type with 
the uniform charge-compensating electron background. Classic and quantum regimes 
of ion motion are considered. The calculations in the classic temperature range agree 
well with previous Monte Carlo simulations. At these temperatures the shear modulus 
is given by the sum of a positive contribution due to the static lattice and a negative 
oc T contribution due to the ion motion. The quantum calculations are performed for 
the first time. The main result is that at low temperatures the contribution to the 
shear modulus due to the ion motion saturates at a constant value, associated with 
zero-point ion vibrations. Such behavior is qualitatively similar to the zero-point ion 
motion contribution to the crystal energy. The quantum effects may be important for 
lighter elements at higher densities, where the ion plasma temperature is not entirely 
negligible compared to the typical Coulomb ion interaction energy. The results of 
numerical calculations are approximated by convenient fitting formulae. They should 
be used for precise neutron star oscillation modelling, a rapidly developing branch of 
stellar seismology. 
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1 INTRODUCTION 



Recent discovery o f quasi-periodic osc i llations (QPO) in soft 



^amma-repeaters lllsrael et al 



20051 : IStrohmaver fc Wattd 



I2OO5I : IWatts fc Strohmaveijl2006 ) may be opening up an ex- 



citing possibility into studying neutron stars by methods 
of seismology. The QPO are thought to be related to neu- 
tron star vibrations and, in particular, originally, they were 
thoug ht to be related to torsion al vibrations of neutron star 
crust (|Duncanlll99i : iPirdbOOsI ). 

Even though it is now understood that the mecha- 
nism of neutron star oscillations is likely more complex 
and involves global oscillations of crust and core, co u pled 
by t he frozen-in magnetic field (|Levinl |2006| . l2007l : lLe3 
|2008| ). it is still possible that the actual oscillation frequen- 
cies are related to pure crustal frequencies, the important 
controlling factor being the mag netic field streng t h and 
geometry jGlampedakis. Samuelsson fc Andersson [200^ : 



The bulk of the neutron star crust is made of fully 
ionized ions (of varying charge Ze and mass M) in crys- 
talline state, immersed in a nearly uniform strongly de- 
generate electron gas. More specifically, the ions form a 
crystal, if the local temperature T falls below the melt- 
ing temperature Tm = Z^e^ /{aVm), where Fm ~ 175, and 
a = {4:Tm/3)~^^^ is the ion sphere radius (n is the ion num- 
ber density, ~ 1). Typically one assumes that the ion 
crystal is of body-centered cubic (bcc) type, as this structure 
is preferable thermodynamically for strictly uniform electron 
background. 

The state of the electron subsystem depends on mat- 
ter density. We shall limit ourselves to such (not too low) 
densities, where electrons are degenerate and ions are com- 
pletely pressure ionized (p > WAZ g cm^ '^, where A is the 



ion mass number; see for discussion Pethick fc Ravenhalll 



IWatts fc Strohmavei 2007 : van Hoven fc Levin 20101 , and 
references therein). The crustal torsional vibration frequen- 
cies are determined by the shear modulus of the solid neu- 
tron star crust. The main purpose of the present paper is to 
calculate this quantity. 

* E-mail:baiko@astro.iofi'e.ru 



ll995l : lHaensel, Potekhin fc YakovlevI 120071 '). At these densi- 
ties the model of uniform charge-compensating background 
of electrons is reasonably good. It gets progressively better 
with the growth of density becoming especially accurate at 
p 3> 10'' g cm"'^, where electrons are ultrarelativistic. The 
crystal of fully ionized ions with the uniform background of 
electrons is known as the Coulomb crystal. 

In inner crust, at densities above the neutron drip den- 
sity pd ~ 4.3 • 10^^ g cm~"^, in addition to the Coulomb 
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crystal of ions and electrons, there are neutrons not bound 
in the atomic nuclei. The details of neutron interactions with 
nuclei are not known very well. It seems plausible that the 
properties of a strictly static crystal, i.e. a crystal with nu- 
clei fixed at their lattice nodes, are determined by Coulomb 
forces. By contrast, the motion of nuclei about the lattice 
nodes may be affected by the presence of neutrons. Lack- 
ing a good model of this effect, we shall assume that it can 
be described by an effective nucleus mass and renormalized 
ion plasma frequency within the framework of the Coulomb 
crystal model. 

The main purpose of this paper is thus to study the 
shear modulus of the Coulomb cr ystal. The g roundwork 
for this problem was laid down by iFuchd l| 19361 ) , who cal- 
culated the shear modulus of the sta t ic bcc Coulomb lat- 
tice. More recently, lOgata fc Ichimarul (|l990l ) calculated the 
shear modulus of the bcc Coulomb crystal taking into ac- 
count the motion of ions about their lattice nodes. In that 
work shear modulus was found n umerically with th e aid 
of Monte Carlo simulations (e.g., iBrush et al.l Il966l ). By 
the nature of the method, the mot ion of ions was treated 
classically. IStrohmaver et al] (Il99ll ) further remarked that 
quantum effects were not important due to the smallness of 
the ion plasma temperature compared to the typical la ttice 
electrostatic energy. Finally, iHorowitz fc Hughtol l|2008l 'l cal- 
culated the shear modulus of the Coulomb crystal taking 
into account weak electron screening in the Thomas-Fermi 
model. This calculation was done numerically using molecu- 
lar dynamics method. Again the motion of ions was strictly 
classic. 

In this paper we shall use the thermodynamic perturba- 
tion theory to find the shear modulus of the Coulomb crys- 
tal with ion motion included in the har monic lattice model 
framework. Unlike numerical methods of lOgata fc Ichimarrj 
l| 19901 ') and iHorowitz fc Hughtol (120081 ) this approach is ca- 
pable of tackling quantum effects. The quantum effects are 
known to be important especially for lighter elements at 
higher densities. Quantum effects in the problem of Coulomb 
crystal elastic moduli are studied for the first time. 

In addition to the outermost envelope of the external 
crust, the Coulomb crystal model fails in the 'nuclear pasta' 
region of the inner crust at densities p > 10^* g cm~^, where 
nuclei become nonspherical. E stimates of the shear modu lus 
in this layer were reported bv lPethick fc PotekhinI (|l998l ). 

Besides neutron star crusts, the Coulomb crystals are 
expected to form in solid cores of white dwarfs, to which the 
present results thus also apply. 



2 GENERAL THEORY 

A Coulomb crystal is composed of ions with charge Ze, ar- 
ranged in a crystal lattice with equilibrium lattice sites Ri, 
immersed in a rigid background of electrons (charge — e). 
Background volume element coordinates are denoted as r. 
Suppose the crystal is deformed uniformly. Then the lat- 
tice remains perfect, but its equilibrium nodes move to new 
locations Xj. The background volume element r shifts to 
a new position x. Using the displacement gradients Ua/3 = 
dXa/dRjj — Saj} (for uniform deformations Uaji — const), one 
can express the coordinates of the new positions via those 



of the old ones as 



(1) 

(2) 



In this case, greek indices denote Cartesian coordinates. We 
do not distinguish between upper and lower indices, and 
summation over repeated greek indices from 1 to 3 is always 
assumed. The potential energy [/ of a uniformly deformed 
Coulomb crystal can be written as 

U 



1 



E 



Xi -f uj — Xj — uj 

n dr n 
+ - 



\Xi + ui — x\ 



dT'idT'2 

\xi - X2\ 



(3) 



where uj is the J-th ion deviation from its deformed equilib- 
rium position Xi due to thermal and zero-point vibrations, 
and n is the ion number density in the nondeformed config- 
uration {R}. Integrations are over the nondeformed crystal 
volume V, prime means that terms with I = J are omitted. 
Since a crystal lattice realizes a local energy minimum with 
respect to small ion deviations from the lattice nodes, the 
energy U can be approximately expressed as 



U 



W\{X}) + SU 



(4) 



where U^^{{X}) is the potential energy of the static de- 
formed lattice (i.e. the energy of the lattice with all ions 
located at the lattice nodes), and SU is the second order 
term of the Taylor expansion in powers of uj: 



SU 



9 ^1.1 



(5) 



dX^idX\ 



5ij — 1 5i.j 

\Xi-Xj\ y^\x,-x. 



n dr 



\Xi 



(6) 



The last term in the curly brackets serves to cancel an in- 
finity arising in the second one. 

The potential energy of the static deformed lattice can 
be expanded in powers of the displacement gradients as fol- 
lows 

:^[/"({X}) = :^(7"*({ii})-|-S'Q;3Ua/3-|-i5'Q/3^AMQ/3U^A. (7) 

In this case, [/^'({i?}) is the potential energy of the static 
nondeformed lattice, while the first order expansion coeffi- 
cient can be expressed via electrostatic crystal pressure: 
S^'a = —P^'^&afi [note, that to first order W jV = Uaa for 
an arbitrary deformation, where SV is th e volume change 
due to the deformation; also see Sect.lSland fWallacel l| 19671 )]. 
Sai3-yX a-re the static lattice elastic coefficients. Elastic coeffi- 
cients as second derivatives with respect to the displacement 
gradients were introduced by Huang (1950). 

Since t he static lattice coefficients are well known 
(|FuchJl936l . see also Sect.©, the main subject of the present 
paper is the temperature and density dependences of the 
elastic coefficients associated with ion vibrations about their 
lattice sites. Thus we shall focus on SU of Eqs. (U) and 
(O. Making use of t he standard in solid-state theory (e.g., 
iBorn fc Huand[l954l ) collective coordinates A'^ (for brevity 
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of notation we consider simple lattices only, i.e. lattices with 
only one ion in the elementary cell): 



1 



VMN 



A'^ exp (ifc • 



Ri 



(8) 



where wavevector k belongs to the first Brillouin zone of the 
nondeformed lattice, while M and A'^ are the ion mass and 
total number, 5U can be written as 



SU 



4N 



2MN 



At fixed J, the sum over / in Eq. ((9|, along with the sum over 
K and the integral in Eq. ((6]), can be extended to infinity. 
The remaining finite sum over J in Eq. ([9]) produces NS^kk/ . 
Then SU can be rewritten as 



SU 



Z e 
2M 



1^0 



1 d'x-^ 



(10) 



where the deformed dynamic matrix D'^"{k, {X}) has been 
introduced. That matrix can also be expanded in powers of 
the displacement gradients and thus related to higher order 
nondeformed lattice energy derivatives: 



D^"{k,{X}) ^ D^"{k,{R}) + D^Jf,{k)u^p 
K^pik) = ^^'[i_exp(ifc.fi)]i?;3 



R 

,3 D-1 



ORcdR^dR^ 



K^xik) = [1- exp (ifc.fi)] J?;3i?, 



R 

d^R-^ 



ORo^dR-ydRudR, 



(11) 



(12) 



(13) 



The prime means that summations over R run over all non- 
deformed lattice vectors except R = 0. Practical expressions 
for the coefficients of the expansion (|lip can be found in 
Appendix A. T hey were obtained us ing the standard Ewald 
technique fe.g.. lBorn fc Huanelll954l '). In Eqs. ((T0|, (O, and 
H13p it is understood that the apparent divergences are can- 
celed by the electron background term of Eq. In practical 
formulae of Appendix A this is refiected by the absence of 
the q — term in sums over reciprocal lattice vectors. 

Nondeformed dynamic matrix D'^" {k, {R\) is given by 
the expression for _D'"'(fe, {X}), stemming from Eq. \1Q\ . 
with Xi replaced by Ri. This matrix determines frequen- 
cies ijjfcj and polarization vectors euj of nondeformed crystal 
oscillations: D^" {k, {R})ekj ~ uj^j^kj^ where j enumerates 
oscillation modes with given fc [j = 1,2,3 for simple lat- 
tices). 

Let us expand over the basis of e^j: A'^ = 

X]f=i6fejQtej- Then the oscillatory potential energy (I10|) 



shall consist of three parts SU 
Ho 



Ho + Hi + H2, where 



Hi 



Ho = 



kj 



-kj' 



kjj' 



-kj' 1 



(14) 
(15) 
(16) 



kjj' 



and, following iBorn fc Huand ||1954| ). we have introduced 
quantities 



E V^j exp [ifc ■ {Ri - Rj) + i(fe + fc') • R.,] . (9) - 



^kj^-kj' 



(fc) = e^je-kj' K'^d-yA'^) 



(17) 

kj^-kj'^afl-,\\"-l ■ (18) 

In this case, Hq is the oscillatory potential energy of the non- 
deformed lattice, while H\ and H2 represent a perturbation 
of this quantity due to the deformation. 

In quantum mechanics coordinates Qkj become opera- 
tors. It is convenient to switch to second quantization rep- 
resentation, where operators Q^j are expressed via phonon 
creation and annihilation operators a^^ and akj'- 



Qkj ~ 



-kj - 



2cjfc 



(19) 



It is now possible to obtain expansion of the free energy in 
powers of the displacement gra dients using the thermod y- 
namic perturbation theory (e.g.. iLandau fc Lifshit j|l980l ): 



= E + E' 



E. 



(0) 



p(0) 



+ 



2r 



+ ■ 



(20) 



In this case, V — Hi + H2 is the perturbation operator, 
indices n and m run over all possible unperturbed quan- 
tum states (which in second quantization means a sum over 
all possible phonon occupation numbers in all modes), and 
Wn = exp{(Fo — En^)/T} is the probability of the quan- 
tum state n, Fo and i?!"' being unperturbed free energy 
and quantum state energy. Terms with n = m in the second 
sum are excluded. 

For simple lattices, e^kj ~ Skj, while all matrices D on 
the r. h. s. of Eq. (|ll|l , are real and symmetric with respect 
to their upper indices fi and v. Consequently, Eq. (|20p can 
be written as {SF)/V = S'^^Uaji + Q-^S^^f^^^^Uaiiu-fX + 
with 



■V ^ • 



h 



2V 



kj 



-T- 



h 



(1 -I- 2nkj 



(21) 



2V ^ 2ujkj 

kj 



1 4- 2nkj) 



+ 2E 



ap 7A 



2V 

kj 



nkj + l)7ife. 



kj' 



(22) 
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where Ukj = [exp (hoj^j/T) — 1] ^ is the average occupation 
number in a Bose system, and the ar gument fc is imp l icit fo r 
all -S's. Note a typo in eq. (41.38) of iBorn fc Huand ljl954 ). 
which differs from our expression H22p by the absence of a 
factor 2 in front of the X^j'^^j square brackets. The term 
with the double sum over j and j' in Eq. (|22p can be also 
written as 



2 

/3 



kj,j'>j 



with (jfej 



2a;, 



-(l + 2nfej) 



which removes the singularity associated with degenerate 
phonon modes ujkj = ^kj' ■ 

The expansion of SF in powers of displacement gradi- 
ents is done at fixed temperature. Accordingly, the elastic 
coefficients of Eq. (|22|l are isothermal. Adiabatic elastic co- 
efficients are defined via expansion of energy in powers of 
displacement gradients at fixed entropy (see Sect. [5}. 

Besides the Huang expansion in powers of the displace- 
ment gradients it is customary to expand thermodynamic 
potentials in powers of the Lagrangian strain parameters 

Vafl = ^{Ua/l + U0a + UXaUXfl) ■ (23) 

In this way one obtains {5F)/V = CafsTjafs + 
0.5Ca/3^Af7Q/3'77A + . . • Since e xpansions in t erms of Uq/j 
and riap must coincide, one has ()Wallacdl 19671 ) 



(24) 



with the same relationships holding for partial contributions 
(e.g., 'st', 'ph', etc). Coefficients C have complete Voigt sym- 
metry, that is Capjx — Ci3a-y\ = C-yXa/s ■ lu general, this 
is not the case for the 5'-coefficients. In cubic symmetry 
there are only three nontrivial C-coefficients dm, C1122, 
and C1212. In Voigt notation they are known as Cii, C12, 
and C44, respectively. 

The results of numerical calculations of elastic coeffi- 
cients are presented in Sect. [S] 



3 RELATION TO PHONON PRESSURE 

As shown, e.g., bv lWallacd (|l967l ). in the presence of an ini- 
tial stress in the nondeformed configuration, —Sap (whether 
isothermal or adiabatic) is equal to this stress. In the case of 
a Coulomb crystal in neutron star crust such initial stress is 
produced by pressure. Consequently, Sap = —PSap, and 
Safi-yx ~ —PSct'ySfix + Cafi-yx- It is thus clear that while 

C1212 = C122I, <S'l212 7^ S1221 ~ <S'l212 + P- 

If the total free energy is a sum of several partial con- 
tributions, their first derivatives yield partial pressures. Just 
like S^/j in Eq. ([7| is related to electrostatic crystal pressure, 
•S'S = -P^'^S^is, where P^^ is the phonon pressure. Phonon 
pressure is found as a volume derivative of the phonon ther- 
modynamic potential Q^^^: 



kj 



+ Tln 



1 — exp 



kj 



T 



(25) 



(26) 



and thermal motion, respectively. The volume dependence is 
contained only in phonon frequencies. The ratios of phonon 
frequencies to the ion plasma frequency Wp = ^ 'InnZ'^e'^ /M 
are universal functions for a given lattice type, and thus 



iUkj (X ujp (X n 



oc y-^/^ It follows, that 



pph 



1 

2V 
1 

4V 



E 

kj 



hjjj, 



kj 



hjjj, 



fcj 



exp [hwkj/T) 



kj 



We can, therefore, assert the following identity: 



E 

hj 



^kj 



+ oJkjSap (1 + 2nkj) = 



(27) 



(28) 



This identity can be proven directly (at least for bcc lattice). 
First, we notice that from Eq. (|17|) together with explicit 
formulae (|A1|) and (|A2|) from Appendix A, it is clear that 
^i:y{ki) = —'l>iy{k2), whcrc fci and fc2 differ from each other 
by the sign of their a;-coordinate (or t/-coordinate) and like- 
wise for ^j^jj with other pairs of indices a ^ p. By contrast, 
^ii, ^"y, and 4>ll do not change under a sign change of any 
of their argument coordinates. If, on the other hand, fci and 
fe2 differ by the interchange of x- and y-coordinates, then 
$ii.{ki) = #^1,(^2), <Pii{ki) = <Pii{k2), and likewise for in- 
terchange of X- and z- or y- and ^-coordinates. This means 
that Ylti ^'Jp{ki) = l6[$iUk)+<Pii{k)+'Pii{k)]S^i3. The 
sum on the 1. h. s. is over 48 Brillouin zone vectors (with 
identical length jfcj), obtained from fc by 6 possible permu- 
tations of absolute values of its Cartesian coordinates and, 
for each permutation, by 8 possible combinations of signs 
assigned to those coordinates. 

For uniform compression SuJkj = —ijJkjSV/{2V) 
and = Sa3SV/{3V). On the other hand, tOkj — 



-kj 



e''_kjD^'' {k,{R}), and therefore, 5uikj = 



kj 

be- 



cause variation of a polarization vector must be orthogonal 
to it in order to maintain its unit length. Combining these 



results we obtain ^>;^^ — 



SiOkj, which proves Eq. 



It is 



The first and second terms in Eq. (|26|) describe zero-point 



obvious from the derivation, that in place of (1 + 2nkj) we 
can have an arbitrary function of | fcj and j. 



4 EFFECTIVE SHEAR MODULUS 

Free energy expansion coefficients, introduced in Sect.[2l also 
determine elastic stress tensor of deformed crystal. If defor- 
mation with displacement gradient Ua/s is applied to a con- 
figuration under initial isotropic pressure P, t he stress tensor 
c"a/3, equal initially to —PS^p, will change by (|Wallacelll967l ) 

SUafi — -Ba/SjxiUjX + UX-i) , (29) 

where 

Ba/B-fX = Sal3^X — P{3axSl3j — ScfiS^x) ■ (30) 

Thus -Biiii — 5*1111, -B1122 = 5ii22 + P, B1212 = S1212 ~ 
B1221 . 

Expression (|29|l allows one to write down linearized elas- 
tic medium equation of motion. In nonuniform matter Saap 
of Eq. (|29p gives Lagrangian variation of the stress tensor. 
In realistic neutron star modelling, equation of motion must 
also take into account magnetic field and non-uniformity of 



© 2011 RAS, MNRAS OOO.nifTOl 



Shear modulus of neutron star crust 5 



matter and ini tial stress, as s ociate d with gravitation [cf., 
e.g., eq. (9) in I Carroll et al.l (|l986l l]. If all these complica- 
ti ons are omitted , the equation of motion reads [eq. (2.23) 
of I Wallace! (|l967l ')] 



pila = B, 



dRadRx 



dRfidRx 



(31) 



where p is the nondeformed mass density and u is the dis- 
placement (so that X = R + u). 
In Fourier space one has 



(32) 



which for cubic symmetry can be expanded as 

22 (— ( /2t2i 2t2i 2t2\ 

pu) u — b\\\\{Uxk^ + Uyky + u^k^) 

+ 2Sii22{uxkxUyky + UxkxUzkz + UzkzUyky) 

+ 2Si22i{uxkxUyky + UxkxUzkz + UzkzUyky) 

+ Si2i2(u'iky + Uyk'^ + u^k'l 

+ u^kl + ulk'i + ulkl) . (33) 

From Eqs. H29p and (|33l) it is clear that in a perfect crys- 
tal it is -B1212 = S'1212 that produces a response to a shear de- 
formation. For this reason we shall call <5'i 2i2 the elastic shear 
coeffic ient. However, it was proposed bv lOgata fc Ichimarul 
l|l99d ) that in order to describe transverse modes in neu- 
tron star crusts, presumably composed of many small ran- 
domly oriented crystalline domains, one has to consider a 
directional average of the above equation. In particular, one 
has to average over directions of u perpendicular to fc and 
then over all possible directions of fc. The resulting isotropic 
phase velocity cj/k should then be equated to effective shear 
wave speed yTfeft/p, where ficn is the effective shear modu- 
lus. The Lagrangian stress tensor variation is then approxi- 
mated by the isotropic medium expression [e.g., eq. (15) in 
ICarroll et all (1986)] with /icff in place of the shear modulus. 

The necessary averaging is easy to carry out. Firstly, 



'all 



kaks 



And secondly, averaging over angles of k yields 



Sii 



5*1221 + 45*12 



(34) 



(35) 



The combination on the r. h. s. of Eq. p5|) is just the effective 
shear modulus fies in question. 

It is not quite clear whether this averaging oversimpli- 
fies the real situation in neutron star crust. Firstly, it is not 
known how small and randomly oriented the crystalline do- 
mains making up the crust really are, and whether or not 
one should consider instea d a mo r e regu lar crystal structure. 
For instance, as shown bv lBaikd (|2009l l, bcc Coulomb crys- 
tal in magnetic field has minimum energy, if it is oriented 
so that the direction of the magnetic field coincides with 
the direction towards one of the nearest neighbors. This ef- 
fect is due to a dependence of zero-point energy on mutual 
orientation of the magnetic field and crystal axes. So, it is 
easy to imagine that during star cooling the crust solidifies 
in such a way that the direction towards a nearest neighbor 
coincides with that of the magnetic field. This will produce 
a large scale ordered crystal structure. Secondly, if we agree 
with the notion of crust as a collection of small randomly 
oriented domains, it is not Eq. (|31|) that has to be averaged, 



but the full equation of motion, which differs from (|3H) by 
the presence of important anisotropics due to magnetic field 
and gravitation. 

Since ^eff contains all elastic coefficients (e.g., Sim 
which is related to bulk compressibility) one may wonder 
whether any other subsystem, besides the ion lattice, con- 
tributes to /icff • If a partial contribution to the energy (or 
the free energy) is a function of particle density only (e.g. 
kinetic energy of the degenerate electron gas), the Huang 
coefficients, associated with it, can be written as 



VS, 



dU 



duafjdu^x 



+ 



duapdu^x dn^ 
drix 9nx d^U 
duafi du-yx drii 



(36) 



In this case, rix is the density in the deformed configuration 
which, to second order in displacement gradients, reads 

: n [1 - TY{Uaf}) 



det(l -I- Uafi) 

+ + U22 + M33 + ^llU22 + ^111^33 

+ U22^i33 + 1tl3lt3l + U12U21 + ^23^*32] , (37) 

where n is the nondeformed density. Consequently, one finds 
F5'iiii = 2n-—+n^-—, 



5*1122 

1^5*1221 
V5*i212 



dU 2 d^U 

dU 
dnx 







(38) 



The important implication is that there are no partial con- 
tributions to neither 5*1212 nor fi^s [see Eq. (|35p ] due to such 
partial contributions to the (free) energy. In particular, nei- 
ther electrons nor d ripped neutrons (in the stan dard model 
of neutron gas, e.g.. IShapiro fc Teukolsk"vlll983t ) in neutron 
star crust contribute to the effective shear modulus. 



5 ISOTHERMAL AND ADIABATIC ELASTIC 
COEFFICIENTS 

In the previous sections we have found formulae for isother- 
mal Huang coefficients 5*^J^^;^ and effective shear modulus 

Moff = (•S'n\i-5fl22-'S'f22i+45*S'2\2)/5- Adiabatic Huang co- 
efficients may be defined in the same way, the only difference 
being that the energy is expanded in powers of displacement 
gradients Ua/s (instead of the free energy). Adiabatic coef- 
ficients are likely much more appropriate for neutron star 
seismology. In this section we show that isothermal and adi- 
abatic 5*1212 as well as /ieff are actually the same. In order 
to prove this, we note that 



dE 



diE,S) d{u^p,T) 



d{Ual3,T) d{Uafl, S) 

dE \ _ (dE\ f dS 
df 



dUa 



dE \ , fdE\ ( dT 
df 



dS_ 

dT 



dUap J J, 



+ 



dUap J g 



(39) 
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Substituting 



fdE_ 



into Eq. (I39|l we see that 
dE 



dE 
'dS 



dS_ 

dT 



(40) 



du. 



Q/3 



dE 
dS 



dE 



+ T 



/ dS 

\dUal3 J J. 
dS 



dE 
du 



a/3 



(41) 



where T = [dE /dS)u^ij was used. (Keeping displacement 
gradients Uaft fixed ensures that the volume and shape of 
the crystal does not change.) Since F — E — TS, Eq. (PTj) 
means that 



dE 

dUc^B 



dF 

dUcB 



(42) 



Therefore 
d^E 



duapdu^x 
d^F 



duapdu^x 



du. 



+ 



( dF ^ 




du^x J 


d 


( dF 


dT 


\duj) 



dT 

dUaH 



(43) 



Expressing derivative of T from Eqs. (|39p and using (|42p we 
obtain 



d^E 



duapdu-fX 



+ T 



s 

dS 



d^F 



duapdu^x 
dS 



dUafI J J, \dUjX 



dE 
df 



(44) 



Since 5* = —{dF/dT)u^f,, and {dF/duap)T oc Sap [i.e. zero 
for a 7^ /3 and same for all a = P; cf proof of Eq. (|28l) ]. we see 
that there is no difference between adiabatic and isothermal 
Huang coefficients 5'i2i2 as well as 51221. Also, we see that 
the difference between adiabatic and isothermal coefficients 
5*1111 is the same as that for coefficients Sii22- This ensures 
that adiabatic and isothermal ^off are the same. 



6 NUMERICAL RESULTS 

In this section we present results of numerical calculations of 
the elastic coefficients for the bcc lattice. For such a lattice, 
one only needs to calculate four coefficients entering Eq. 
((35)) . 5iiii, S1122, 5i22i, and 5i2i2. AH the other coefficients 
with 2 pairs of identical indices are equal to these ones (e.g., 

5llll — 52222 = 53333, 5l212 ~ 52121 = 5l3l3 — . . . , 5l221 = 

52112 = 5i33i = . . . , 5ii22 = 52211 = 5ii33 = . . .), while 
the rest of the coefficients are zero. Alternatively, one can 
calculate 5iiii, 5ii22, 5i2i2, and pressure P. In what follows 
we shall focus on static lattice and phonon contributions. As 
discussed in Sect. [J] there are no other contributions to the 
shear modulus. 

The static lattice elastic coefficients are well-studied in 
the literature. Practical expressions for the coefficients of 
the expansion Q can be found in Appendix B. They are 
derived using standard Ewald technique and given here for 
completeness. The numerical results obtained using these 
expressions ar e given in T able [1] These values agree with 
the results of iFuch j (|l936l ). who calculated lattice energy 



expansion coefficients for two types of elastic deformations, 
A = -P=* + 5fiii - 5fi22 and 2B = 31^2, for bcc and 
face-centered cubic lattices (see also I Wallace! 1967h . 



Computation of phonon coefficients H22p requires inte- 
gration over the first Brillouin zone, 

Methods of such integration have been developed elsewher e 
(|Albers fc Gubernatid[l98ll : lBaikdl2000l : iBaiko et aLlboOll l. 

It is sufficient to integrate only over the primitive cell of 
the Brillouin zone, described, for bcc lattice, by inequalities 
kx ky ^ kz ^ and kx + ky ^ 2Ti/ai, where ai is the bcc 
lattice constant: naf = 2. The phonon frequencies are even 
functions of k^, ky, kz and are invariant under an arbitrary 
permutation of kx, ky, k^. In order to formulate symmetry 
properties of other quantities entering (I22|l we denote an 
arbitrary permutation of {x,y,z) as {p,a,T). If one inter- 
changes k„ ^ kr then ^;^^pp, ^H^t (equal to ^i^rra 
(equal to ^5;^^^^) remain the same, whereas 'P-'Jaaa ^ ^tttt, 

^ ppaa ^ ^ppTT'i ^paap ^ ^pTTp: ^ pa per ^ ^ pT pr ^ ^apap ^ 

'^t'ptpj and <^CTro-T ^ ^TCTTCT- Thc same relationships gov- 
ern symmetry properties of products i^J^ , viewed as 
4-index {ajSjX) quantities. Additionally, ^H-tt = ^H-ra, 
and = ipii, so that <?ji#^i = ^H^H = il>i.{^li. 
Both '^''^fj^x '^^'^ products (l>"p <l>ij'^ are even functions of 
kx , ky , kz ■ 

Let us give the recipe to calculate an arbitrary 5^*^ co- 
efficient from Eq. (|22|l . In order to obtain, for instance, Sf^n 
one has to integrate over the primitive cell of the Brillouin 
zone the integrand in Eq. (|22|) with a/3-yX — xxxx, yyyy, 
and zzzz. Then add the three results together and multiply 
by 16 to account for the remaining two permutations for 
positive kx, ky, kz and for 8 possible combinations of signs 
of kx, ky, kz- The same recipe is applied in the case of 55'J'22 
and 5^^221 • 



The situation with 5^212 is more complex because, for 
instance, ^lyxy 7^ ^l^yx- This requires integration over the 
primitive cell of the integrand of Eq. (|22|l with afi-yX — 
xyxy, yxyx, xzxz, zxzx, yzyz, and zyzy, addition of all 6 
of them together and multiplication by 8. 



We shall now describe the results of numerical calcula- 
tions. In Fig. [T]we show —5^212 (upper dashed curve) and 
"Meff (lower dashed curve) in units of nTp as functions of 
T/Tp, where Tp = hUp is the ion plasma temperature. Quan- 
tities 5^212 and /i^^ are negative. Thus they reduce the re- 
spective static lattice values (Table [l| and weaken lattice 
resistance to the shear strain. Dots show the same quanti- 
ties with quantum effects explicitly excluded. In this case 
the elastic coefficients are always proportional to T. These 
results are obtained by setting fikj = T/hio^j and retaining 
only the highest order terms in T in Eq. ([22}. At higher 
temperatures the classic curves merge with the exact re- 
sults, whereas at lower temperatures the quantum effects 
dominate. The exact curves do not decrease beyond certain 
values corresponding to ion zero-point motion. In this tem- 
perature regime the perturbation theory has clear advantage 
over Monte Carlo or molecular dynamics methods, in which 
the ion motion is treated classically. 
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Table 1. Static lattice clastic coefficients (bee) in units of nZ^e^/(2a; ), where a; is the bee lattice constant: naf 



c-st 

'-'1111 


'-'1122 


est 
"^1212 


^st 


A 




-1.4848079 


-0.47067387 


0.74240395 


1.2130778 


0.19894377 


0.48523113 



10 



H 

n 

3. 
I 
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0.1 



1 1 1 — I — I I I I I 1 1 1 — I — I — r-n 




classic 



_L 



10 



10 
T / T, 



Figure 1. Phonon elastic coefficient —3^212 (upper dashed curve) 
and effective shear modulus —fJ^^g (lower dashed curve) in units 
of nTp vs. T/Tp. Solid curves and dots show analytic fits [Eqs. 
1145 I I and II46| |] and classic numerical results, respectively. 
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Figure 2. Thermal contributions — Mgff (solid line) and —S\^^2 
(dashed line) and zero-point contribution to — (dots) in units 
of nTp vs. T/Tp. 



We were able to fit the phonon shear coefficients as 

1/3 



ph 



~nTp 



0.3686^ + 136.6 ( ^ 



0.5903^ + 439 ,' J 



(45) 



1/3 



with the dashed ones). They reproduce exactly the classic 



,ph 



fit is 



and zero-point limits. The maximum error of the /i^g 
2.2% at T/Tp ^ 0.06. The maximum error of the 8^212 fit is 
1.0% at T/Tp ^ 0.08. 

If one subtracts the T -> limit (i.e. the zero-point con- 
tribution) from the phonon elastic coefficients, the remain- 
ing part will correspond to thermal ion motion (we denote 
this part by a superscript 'th'). In analogy with the deriva- 
tion of the Debye T'^-law for specific heat, one can show 
that at low T this thermal contribution to the elastic co- 
efficients behaves as (e.g., Bor n fc Huan g 1954 ). In the 
classic regime of high T, both 51212 and fi^ are, naturally, 
proportional to T (cf Fig. [T]). Our numerical calculations of 
— and —5*12 reproduce both asymptotes and are shown 
in Fig. [2] by solid and dashed lines, respectively. Dots show 
the zero-point contribution to /i^^. Since it is the sum of 
the zero-point and thermal contributions that make up the 
total phonon elastic coefficient, the part of the thermal 



Table 2. Asymptote parameters /, g, and h in units of nTp, where 
/(T/Tp)'* is the thermal quantum asymptote, gT/Tp is the classic 
asymptote, and h is the T = (zero-point) value of — and 



-,ph 
'l212- 







/ 


9 


h 


(46) 


ph 


1.43 X 10-* 


5.14 


0.369 




oph 
■^1212 


1.35 X 10-* 


7.60 


0.590 



contribution is practically always negligible. We summarize 
parameters of various asymptotes in Table [21 

Finally, let us compare our results with those of other 
authors. In Fig. [3] we show various approximations to the 
total effective shear modulus in units of nZ^e^/a versus 
r = Z^e'^/(aT) for fully ionized ^^C at density 10^° g cm'^ 
In this case a — (47rn/3)~^'''' is the ion-sphere radius and 
r is the standard Coulomb coupling par ameter. Bars repre- 
sent original Monte Carlo calculations of lOgata fc Ichimarul 
(199( |), while dash-dott e d cur ve shows the fit to these data 
from IStrohmaver et al.l (|l99ll ). The solid line is /i^g -I- /i^^, 
where (1% is from Table[T]and fi^^ is given by Eq. (|45|) . Dots 
show the sum of and the classic asymptote of /i^^ (dots 
in Fig. [T|. This curve thus represents results of a purely 
classic calculation and may be directly compared with the 
Monte Carlo study. The dashed line shows results of molec- 
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Figure 3. Total effective siiear modulus in units of nZ^e^ /a vs. 
r for fully ionized ^^C at density 10^" g cm~^. Solid line and dots 
represent results of our quantum and classic calculations, respec- 
tively, bars show Monte Carlo (MC) data, dash-dotted curve is 
the fit to them, and dashed curve represents molecular dynamics 
(MD) results. 



ular d ynamics simulations reported bv lHorowitz fc Hughtol 
l|2008[) . 

First of all, we note that our class i c calc ulation (dots) 
matches the bars of lOgata fc Ichimarul l|l990t ) with the ex- 
ception of the one at F = 200, where there is a discrepancy 
of about 3% between the dotted curve and the upper tip 
of the error bar. One possible reason for this discrepancy is 
associated with anharmonic corrections to the free energy, 
not taken into account in our perturbative calculation. At 
F = 200 the ratio of anharmonic and harmonic energies 
of the classi c crystal \s (AiNT /T) / {2,NT) « 1.8%, where 
Ai » 10.84 (|P ubin"l990|) . On t he other hand, it is not clear 
from lOgata fc Ichimaru l|l990l ) what is the confidence level 
of their error bars. 

The difference between the quantum calculations (solid 
curve) and the classic ones depends on temperature and 
reflects the results shown in Fig. [1] The deviation of the 
quantum curve is mostly due to zero-point ion vibrations 
and is the strongest at lower temperatures (higher F's). 
The relative magnitude of the deviation is proportional 
to ahwp/{Z'^e^) ^ 0.05pi^V(^6Aif ), where Ze = Z/6, 
Ai2 = A/12, and pio is the mass density in units of 10^" 
g cm"'^. For carbon at p — 10^'' g cm~^ and low T's, the 
classic efi'ective shear modulus exceeds the more accurate 
quantum result by up to 18%. 

As expected, at higher temperatures the quantum curve 
merges with the classic one. It is well known that for a 
Coulomb system the liquid state is energetically preferable 
to the crystal at F < 175. If there is a crystal at these 
F's, it is in the metastable overheated state. We have ex- 
tended our calculations into this temperature range for the 
sake of a qualitative discussion. Obviously, the further we go 
into the overheated crystal regime, the less accurate our re- 



sults become due to the growing importance of anharmonic 
effects. In accordance with the general picture of Fig. [TJ 
where the phonon contribution into the effective shear mod- 
ulus is negative and grows by absolute value with tempera- 
ture, the total effective shear modulus drops and, eventually, 
reaches zero at F ~ 45. This point should be regarded as a 
lower bound for the disappearance of the shear modulus in 
a Coulomb system. The actual crossing point is expected to 
occur at higher F's and will be sensitive to the anharmonic 
effects. 

Molecular dynamics calculations of the effective shear 
modulus in a Coulomb system with compressible elec- 
tron background, produci ng screening, were performed by 
iHorowitz fc Hughtol (|2008l ). These results are shown in Fig. 
[3] by the dashed line. The electron screening was described 
by the Thomas-Fermi (TF) model. By its nature the molec- 
ular dynamics approach does not take into account quantum 
effects. We observe a decrease of the shear modulus (with 
respect to the classic dotted curve) reflecting the fact that 
the crystal with screened Coulomb potential is less tightly 
bound than the pure Coulomb crystal. 

It is well known that to lowest order in electromagnetic 
coupling, the screening by relativistic degenerate electrons 
is des cribed by the (RPA) dielectric function of I Jancovicil 
(|1962| ). The resulting effective ion-ion potential is more com- 
plex than the simple Yukawa potential arising in the TF 
model. Electron screening affects both properties of the 
static lattice and the phonon modes. As shown by iBaikol 
(2002), the TF model overestimates the correction to the 
lattice electrostatic energy as compared to the more accu- 
rate RPA model. By contrast, the properties of the phonon 
modes obtained in the TF and RPA models are very similar. 
Since the main contribution to the effective shear modulus 
comes fro m the electrostatic energy, we expect, that calcu- 
lations of iHorowitz fc Hughtol (|2008l ) in the TF model over- 
estimate the importance of the electron background com- 
pressibility. It would be interesting (and is not difficult) to 
calculate the electrostatic contribution to the effective shear 
modulus using the full RPA model. 



7 CONCLUSIONS 

We have calculated elastic shear coefficient 51212 and ef- 
fective direction averaged shear modulus ficS for neutron 
star crust matter. Both coefficients are sums of a static lat- 
tice term and a term originating from the motion of ions 
about their lattice nodes. While the static lattice contribu- 
tions (Table [l| were well known previously, only numerical 
simulations existed for the ion motion terms. Using ther- 
modynamic perturbation theory, we have expressed the ion 
motion terms via integrals over the first Brillouin zone of 
quantities given by rapidly convergent lattice sums, Eq. (|22|) . 
The integrals were then evaluated numerically and the re- 
sults were fitted by simple analytic formulae, Eqs. (|45|l and 
(gSll. 

The main advantage of the numerical methods (Monte 
Carlo and molecular dynamics) is their ability to include 
anharmonic effects to all orders. The main advantage of 
the perturbation theory is its ability to include quantum ef- 
fects (within the framework of the harmonic lattice model), 
greater transparency of the results and much lesser computer 
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time requirements. The anharmonic effects can also be taken 
into account in the perturbative approach, but that would 
make all the equations much more cumbersome even if only 
the lowest order anharmonic term is retained. Summarizing, 
we can say that numerical simulations and perturbation the- 
ory ideally complement each other as well as serve for mutual 
verification. 

If quantum effects are included, one finds that the ion 
motion contribution can be decomposed into two parts. One 
of them corresponds to zero-point ion motion and is inde- 
pendent of T. The other corresponds to thermal ion motion 
and is cx T* at low temperatures and cx T at high temper- 
atures. The high T asymptote is what one obtains, if the 
calculation is purely classic. At low temperatures the ther- 
mal term is negligible compared to the zero-point term, and 
thus the r* asymptote seems rather unimportant. Our fit- 
ting formulae 1)45^ and (|46p reproduce exactly the high T 
and zero-point limits. 

If quantum effects are excluded, our results agree well 
(c f. dots and bars in Fig. IJJ with Monte Carlo simulations 
of lOgata fc Ichimarul (|l990D . The only discrepancy of about 
3% occurs at r = 200, where anharmonic effects are at their 
strongest (and also the error bars of the Monte Carlo simu- 
lation are at their largest). If quantum effects are included, 
then the main difference with the Monte Carlo results is 
due to the zero-point contribution to the ion motion term. 
Compared to the total shear modulus, that also includes the 
static lattice part, this contribution is important for lighter 
elements at higher densities, where the ion plasma temper- 
ature is not entirely negligible with respect to the typical 
Coulomb ion interaction energy. 

We have demonstrated that neither 51212 nor /icff have 
any contributions from subsystems whose partial free ener- 
gies are functions of particle number density only; that both 
coefficients are the same whether they are evaluated at con- 
stant temperature or entropy; we have also proven an iden- 
tity linking phonon pressure with a coefficient of dynamic 
matrix expansion in powers of displacement gradients, Eq. 

The results, reported in this paper, also apply to crys- 
tallized matter in white dwarf cores. 
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APPENDIX A: 

Both _D^J^(fc) and D^'^^^^{k) are sums of two series, over di- 
rect and reciprocal lattice vectors (denoted R. and G respec- 
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tively); D'";. = D^"™ + D^"'^', D"" , = D^"'-"' + D''"'^' 

' Q/3 a/5 Q/3 ' ctp'^X ctp^X apyX 
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SRaRp Ri_tRv 
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16 
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where 

^(2)g(2) ^ 5aA (57/39^9^ + ^m/S'Jt'?!^ + ^>'/397I7m) 

+ [Sapq^qu + 5nf3qaqv + S^pqaq^i) 

+ 5i.A {Sapq^q^i + Sy^Qaq^ + 5^pqaq-i) , (A5) 

(Sg'*^ = Saxqpq^qfj.qi' + S^xqaq^qt^qu + S^xqaq/sq^iqv 



J6) 



+ S^xqaqpq-iqu + Svxqaqpq-iqii 

+ Sapq-^q^^q^qx + Sypqaqf^qvqx 

+ 5^ii3qaq~,qvqx+ 5upqaq-,qtiqx , 

= qaqpq-^qxqiiqv ■ 



(A6) 
(A7) 



In the above equations, as well as in Appendix B, erfc is 
the complementary error lunction, and A is an arbitrary 
parameter with units of inverse length chosen so that sums 
over direct and reciprocal lattice vectors are equally rapidly 
convergent. Integrals over dq in Eqs. (|A2p and (|A4|) are kept 
for brevity of the formulae. They can be easily evaluated 
with the aid of the 5-functions. 



APPENDIX B: 

Similarly to Appendix A, S^g = S^^' + sf^ and S'^j^;^ 
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